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In  the  present  work  a  physically  based  model  of  direct  methanol  fuel  cell  anode  impedance  has  been 
developed  and  validated  at  different  operating  current  densities.  The  proposed  model  includes  the  two- 
phase  mass  transport  of  both  methanol  and  water  through  diffusion  and  catalyst  layers  and  the  methanol 
oxidation  reaction  involving  CO  adsorbed  intermediate.  Model  simulations  are  in  good  quantitative 
agreement  with  experimental  observations  and  permit  to  evaluate  the  origin  of  anode  impedance  fea¬ 
tures.  Model  results  confirm  that  the  high  frequency  45°  linear  branch  is  caused  by  proton  transport 
limitations  within  the  catalyst  layer  and  that  the  low  frequency  inductive  behavior  is  due  to  surface 
coverage  by  CO  reaction  intermediate.  Moreover  model  predictions  elucidate  the  contribution  to  the 
impedance  of  mass  transport  phenomena  through  diffusion  layer,  that  is  relevant  even  at  low  current 
density  and  increases  along  the  channel  length.  In  particular  liquid  convective  fluxes  are  considered  as  a 
process  of  pressure  buildup  and  breakthrough  at  diffusion  layer  intersecting  fibers,  resulting  in  a 
discontinuous  phenomenon.  By  means  of  this  intermittent  description  it  is  possible  to  correctly  repro¬ 
duce  mass  transport  limitations  through  diffusion  layers,  that  manifest  themselves  as  a  second  arch 
superimposed  to  the  first  one,  peculiar  of  kinetic  losses. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  Electrochemical  Impedance  Spectroscopy  (EIS)  is  a  powerful 
in-situ  measurement  technique  that  permits  to  evaluate  the  kinetic 
and  transport  phenomena  of  an  electrochemical  system  [1,2],  It 
consists  in  perturbing  the  fuel  cell  operation  with  a  small  AC  cur¬ 
rent  signal  over  a  wide  range  of  frequencies  and  in  measuring  the 
voltage  response.  The  module  and  the  phase  shift  between  voltage 
response  and  current  perturbation  are  due  to  the  electrical 
impedance  of  the  fuel  cell  and  provide  useful  information  on  the 
entity  of  internal  fuel  cell  losses. 

Despite  the  potentialities  of  this  measurement  technique  the 
interpretation  of  experimental  observations  is  very  complicated 
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and  modeling  plays  a  fundamental  role  in  the  analysis  of  experi¬ 
mental  data.  Until  now  the  interpretation  has  mostly  been  carried 
out  by  means  of  equivalent  circuit  method  (ECM)  [3—6],  Even 
though  simple  and  fast,  this  method  is  not  reliable,  since  the 
equivalent  circuit  is  not  unique;  moreover  ECM  provides  only  few 
useful  qualitative  information. 

M.  Orazem  [7]  proposed  an  innovative  approach  for  a  quanti¬ 
tative  interpretation  of  EIS  measurements:  the  developed  meth¬ 
odology  consists  in  expressing  the  equivalent  circuit  element  as  a 
function  of  the  physical  parameters  of  the  system.  In  Ref.  [8]  this 
method  was  used  for  the  interpretation  of  cathode  impedance  of  a 
polymer  electrolyte  membrane  fuel  cell  (PEMFC),  but  in  principle  it 
could  be  applied  to  any  electrochemical  system,  including  direct 
methanol  fuel  cell  (DMFC).  However  the  results  still  depend  on  the 
considered  equivalent  electric  circuit.  A  similar  methodology  has 
been  adopted  by  T.S.  Zhao  et  al.  [9]  to  investigate  the  simultaneous 
oxygen  reduction  and  methanol  oxidation  at  the  cathode  of  a 
DMFC. 
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An  alternative  approach  is  represented  by  physically  based  EIS 
models,  that  in  the  last  years  attracted  increasing  scientific  interest 
[10—19].  The  development  and  experimental  validation  of  this  type 
of  models  are  very  complicated,  but  the  results  are  not  related  to 
the  choice  of  a  suitable  equivalent  electric  circuit  and  model  pre¬ 
dictions  permit  to  investigate  the  origin  of  different  impedance 
features. 

In  the  literature  many  works  regarding  both  solid  oxide  fuel 
cells  (SOFCs)  [10—12]  and  PEMFCs  [13,14]  technology  can  be  found. 
Bessler  et  al.  carried  out  an  extensive  physically  based  model  ac¬ 
tivity  on  SOFC  [15—17]:  the  flexibility  and  modularity  of  the 
developed  models,  including  detailed  elementary  kinetic  electro¬ 
chemistry  and  diffusion  processes,  allowed  the  assignment  of  the 
origin  of  SOFC  impedance  features  with  high  accuracy.  Springer 
et  al.  [18]  developed  one  of  the  first  numerical  models  for  PEMFC 
cathode  impedance:  the  first  impedance  loop  was  attributed  to  the 
effective  charge-transfer  resistance  and  double-layer  charging, 
while  the  second  one  was  related  to  the  mass-transport  limitation 
in  the  gas  phase.  Guo  and  White  [19]  presented  a  PEMFC  cathode 
impedance  model  consisting  of  many  flooded  spherical  agglom¬ 
erate  particles,  permitting  a  more  detailed  comprehension  on 
impedance  behavior. 

In  DMFCs  the  anode  catalyst  layer  (CL)  internal  losses  are  not 
negligible,  as  in  PEMFC  technology;  moreover  the  mass  transport 
phenomena  at  anode  side  are  very  complex,  due  to  the  two-phase 
and  multi-component  nature  of  the  flows.  Thus  the  interpretation 
of  impedance  behavior  is  even  more  complicated  than  in  PEFC  and 
SOFC  technology,  in  fact  in  the  literature  only  few  modeling  ana¬ 
lyses  can  be  found  [20—24],  generally  concerning  the  cathode. 

Kulikovsky  reported  a  clear  and  detailed  explanation  of  the 
mathematical  development  of  DMFC  cathode  impedance  model 
[20,21  ],  analyzing  the  effects  of  physical  phenomena  on  the  shape 
of  cathode  impedance  spectrum.  However  the  models  are  not  in¬ 
tegrated  along  channel  length  and  no  attempts  to  fit  experimental 
data  have  been  done.  Sundmacher  et  al.  presented  a  detailed 
analysis  of  the  methanol  oxidation  kinetics  on  a  DMFC  anode  in  a 
cyclone  flow  cell  [23],  but  the  homogeneous  concentration  distri¬ 
bution  over  the  membrane  electrode  assemblies  (MEA)  does  not 
permit  to  investigate  mass  transport  limitations.  Therefore  in  the 
literature  limited  effort  has  been  dedicated  to  develop  a  whole 
DMFC  anode  impedance  model  and  a  detailed  comprehension  of 
the  physicochemical  phenomena  regulating  anode  operation  is 
required  to  further  optimize  components,  operation  strategies  and 
lifetime. 

This  work  proposes  a  complete  and  validated  physical  model  of 
DMFC  anode  impedance,  increasing  the  understanding  of  anode 
operation  and  providing  a  quantitative  interpretation  of  the 
experimental  measurements.  The  work  is  organized  as  follows.  In 
Section  2  the  experimental  measurement  of  DMFC  anode  imped¬ 
ance  is  reported,  subsequently,  in  Section  3,  the  model  develop¬ 
ment  is  described.  Then,  in  Section  4,  the  model  results  and  the 
origin  of  different  impedance  features  are  illustrated  and  finally 
some  conclusions  are  given  in  Section  5. 

2.  Experimental 

2.1.  Cell  hardware 

The  experimental  measurements  were  conducted  using  the 
same  fuel  cell,  as  well  as  equipment,  which  has  already  been 
characterized  in  terms  of  anode  polarization,  performance,  meth¬ 
anol  cross-over  and  water  transport  in  previous  publications  by  the 
authors  [25—27]. 

The  MEA  was  purchased  already  assembled  by  balticFuelCells 
GmbH  with  an  active  area  of  22.1  cnT2,  consisting  of  Pt— Ru  anode 


(4  mg  cm-2,  Pt:Ru  =  2),  a  Nation®  117  membrane  and  Pt  cathode 
(4  mg  cm  2).  Anode  diffusion  layer  presents  gas  diffusion  layer 
(GDL)  without  micro  porous  layer  (MPL),  while  the  cathode  one  has 
got  MPL:  for  this  reason  this  MEA  is  named  MEA  GM,  as  in  previous 
publication  [27],  The  fuel  cell  is  thermostated  at  333  K,  the  anode  is 
supplied  with  a  3.25wt.%  aqueous  solution  of  methanol  at  1  g  min-1 
flow  rate,  while  the  cathode  operates  with  a  hydrogen  flow  rate 
(see  paragraph  2.2)  of  2.14-10-4  g  min-1. 

2.2.  Electrochemical  measurements 

In  hydrogen— air  fuel  cells  contributions  of  the  anode  are  usually 
negligible  due  to  fast  kinetics  of  hydrogen  oxidation  reaction  and  a 
standard  practice  consists  in  neglecting  anode  impedance  [18],  In 
DMFC  the  slow  methanol  electro-oxidation  does  not  permit  to 
clearly  distinguish  anode  and  cathode  contributions  by  measuring 
the  full  fuel  cell  impedances.  However  it  is  possible  to  eliminate 
contributions  of  the  cathode  in  a  half-cell  DMFC  by  feeding  the 
cathode  with  hydrogen,  so  that  protons  are  reduced  and  hydrogen 
is  evolved  [28—31],  In  this  configuration  the  cathode  works  as  a 
dynamic  hydrogen  electrode  (DHE)  and  it  is  suitable  as  a  reference 
and  counter  electrode  for  DMFC  anode  impedance  measurements. 

The  anode  spectra  were  recorded  with  an  Autolab  PGSTAT  30® 
provided  with  a  frequency  response  analysis  module.  The  imped¬ 
ances  are  measured  under  galvanostatic  control  and  the  amplitude 
of  the  sinusoidal  current  signal  is  adjusted  so  that  the  potential 
amplitude  does  not  exceed  10  mV.  The  frequency  is  included  be¬ 
tween  10  kHz  and  50  mHz  with  a  logarithmic  distribution.  The 
obtained  experimental  values  are  processed  by  a  retrospective  use 
of  Kramers— Kronig  transforms  [1,32  [  in  order  to  verify  the  validity 
of  the  measurements.  The  impedance  values  that  do  not  satisfy 
such  relations  are  not  considered  meaningful.  Fig.  1  reports  two 


Fig.  1.  Anode  Nyquist  plots  at  different  current  densities:  a)  raw  data;  b)  after  the  use 
of  Kramers-Kronig  relations. 
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anode  spectra  measured  at  0.075  A  cm~2  and  0.15  A  citT2  before 
and  after  the  use  of  Kramers— Kronig  relations.  It  is  evident  that  the 
high  frequencies  region  is  not  consistent  with  these  relations, 
because  of  instrument  limitation,  and  therefore  the  corresponding 
measurement  points  are  eliminated  from  the  spectra. 

The  anode  impedances  of  Fig.  1  are  an  elongated  semicircle:  this 
is  partially  due  to  proton  transport  losses  in  the  CL,  that  manifest 
themselves  as  almost  45°  linear  branch  at  high  frequencies.  By 
varying  the  operating  current  density,  two  relevant  differences  are 
appreciable  in  the  low  frequency  region.  At  0.075  A  citT2  an 
inductive  behavior  is  evident:  this  could  be  due  to  the  catalyst 
surface  coverage  by  adsorbed  reaction  intermediates.  At 
0.15  A  citT2  the  inductive  behavior  is  no  longer  present,  but  mass 
transport  limitations  seem  to  appear  in  the  low  frequencies  region. 
This  impedance  feature  is  coherent  with  the  increased  mass 
transport  fluxes  through  the  GDL. 

However  in  the  literature  the  interpretation  of  DMFC  anode 
impedance  is  not  fully  consolidated  [33]  and  in  the  next  sections  a 
physically  based  impedance  model  of  DMFC  anode  is  presented  and 
discussed,  in  order  to  provide  a  more  solid  and  quantitative  inter¬ 
pretation  of  experimental  observations. 

3.  Model  development 


still  regulated  by  water  transport  through  the  membrane  and 
cathode  diffusion  layer,  coherently  with  the  previously  pro¬ 
posed  interpretations  of  water  management  [27]; 

•  the  contribution  of  methanol  liquid  convection  through  anode 
GDL  is  considered; 

•  liquid  saturation  is  no  longer  constant  across  the  anode  diffu¬ 
sion  layer  thickness  due  to  the  liquid  convective  fluxes. 


The  equations  that  fully  describe  mass  transport  phenomena 
through  anode  GDL  are  the  following: 


“dgdl,ch3ohsgdl- 


•a /'-•LiUL.L 

OLCH3OH  nG  M  «.  \  OLCH3OH 
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dx 
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In  this  paragraph  a  previously  developed  1D+1D  DMFC  model 

[27,34]  has  been  firstly  integrated  with  a  detailed  description  of  9^gdl°H  ^ChJoh  -  ^ctpoH 

mass  transport  phenomena  through  anode  GDL  and  a  CL  with  a  ax  —  — £gdl  sgdl  ^  £gdl  (  —  sgdl)  ^ 


finite  thickness.  Subsequently,  considering  the  same  mathematical 
approach  presented  in  Refs.  [20,21],  a  DMFC  anode  impedance 
model  has  been  developed. 

3.1.  Diffusion  layer  model 


(3) 
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In  the  ID  +  ID  DMFC  model  presented  in  Refs.  [27,34],  due  to 
the  high  availability  of  water  at  the  anode,  the  phenomena  gov¬ 
erning  water  transport  through  anode  GDL  are  omitted;  the  water 
transport  through  the  MEA  is  assumed  to  be  regulated  by  water 
transport  through  the  membrane  and  cathode  diffusion  layer. 
Moreover  in  the  original  model  the  methanol  flux  through  the 
anode  GDL  is  governed  by  gas  and  liquid  diffusion,  assuming  that 
liquid  saturation  is  constant  across  the  diffusion  layer  thickness  (Eq. 
(18)  in  Ref.  [34]). 

In  the  developed  model  the  following  improvements  have  been 
introduced: 


•  the  water  is  transported  through  the  anode  GDL  by  liquid  con¬ 
vection,  Fig.  2:  however  the  water  transport  through  the  MEA  is 


Fig.  2.  Water  fluxes  through  the  DMFC  during  anode  polarization. 


Eq.  (1)  is  the  Fields  law  of  diffusion  for  two-phase  methanol,  in 
which  the  second  term  in  the  right  hand  side  of  the  equation  is  the 
contribution  of  methanol  liquid  convection.  As  proposed  in  Refs. 
[34],  the  methanol  concentrations  in  gas  and  liquid  phases  are 
considered  in  equilibrium,  described  by  Henry’s  law: 


rGDL,L 

lch3oh 


r/  /-GDL,G 

Kh.ch3ohLCHjOH 


■RT 


(5) 


Eq.  (2)  is  the  governing  equation  of  water  transport,  that  is 
described  by  the  Leverett  function  [35].  Eqs.  (3)  and  (4)  express  the 
mass  balance  of  methanol  and  water,  respectively;  the  liquid  phase 
accumulation  term  is  proportional  to  the  liquid  saturation,  while 
the  gas  phase  one  is  proportional  to  the  void  fraction. 

Particular  attention  has  been  given  to  the  understanding  of  the 
mechanisms  regulating  liquid  convective  fluxes.  In  fact  the  GDL  is 
not  assumed  to  be  always  flooded  with  fully  liquid  pathways,  as 
widely  accepted  in  the  literature  [35,36].  Recent  studies  on  liquid 
convection  through  GDLs  [37,38]  showed  that  liquid  paths  are 
intermittent  and  breakthrough  locations  change  with  time.  These 
flow  visualizations  [37,38]  suggest  that  the  liquid  transport  within 
the  GDL  is  a  process  of  capillary  pressure  buildup  and  breakthrough 
at  the  interface  of  GDL  intersecting  fibers,  Fig.  3  a.  When  the  liquid 
pressure  exceeds  the  breakthrough  value,  the  liquid  expands 
rapidly  through  this  interface  until  the  fluid  contacts  the  next  fiber 
intersections  interface,  Fig.  3b.  Then  the  process  of  pressure 
buildup  and  breakthrough  begins  again  until  the  fluid  reaches  the 
GDL-CL  interface,  Fig.  3c.  At  this  interface,  when  the  capillary 
pressure  exceeds  the  threshold  value,  the  liquid  is  quickly  removed 
from  GDL  pores,  that  become  empty,  Fig.  3d.  The  dynamic  of  liquid 
emergence  and  detachment  from  GDL  surface  is  very  fast,  in  the 
order  of  few  milliseconds,  and  it  is  enhanced  by  GDL  hydropho- 
bicity  [39,40], 
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The  proposed  description  of  intermittent  liquid  convective  flux 
implies  that  the  effect  of  liquid  accumulation  terms  in  the  GDL,  Eqs. 
(3)  and  (4),  is  only  present  meanwhile  liquid  flux  to  the  CL  effec¬ 
tively  occurs,1  or  rather  for  frequencies  higher  than  tens  of  Hertz.  In 
this  work  a  threshold  value  of  75  Hz  is  assumed. 


r 


a7co 

pt'^r 


i  .  i  . 
~  4 T'1  2 T'2 


where: 


(11) 


3.2.  Catalyst  layer  model 

In  the  original  1D+1D  DMFC  model  [27,34]  the  CL  is  assumed 
to  be  monodimensional  in  the  direction  of  channel  length.  In  this 
work  a  CL  with  a  finite  thickness  is  introduced  in  the  model  in 
order  to  take  into  account  proton  transport  limitations.  Moreover, 
coherently  with  the  presence  of  low  frequency  inductive  loop 
(Fig.  1 ),  a  reaction  mechanism  involving  intermediate  adsorbates  is 
considered.  In  particular  the  following  two-steps  reaction  mech¬ 
anisms  with  a  single  adsorbed  intermediate  is  proposed 
[29,41,42]: 

CH3OH  +  Pt  - ^  Pt  -  CO  +  4H+  +  4e-  , 

(6) 

Pt  -  CO  +  H20  - ^  Pt  -  C02  +  2H+  +  20- 

In  the  first  reaction  step  the  adsorbed  intermediate  CO  is 
formed  in  the  process  of  methanol  dehydrogenation,  while  in  the 
second  reaction  step  adsorbed  CO  undergoes  further  oxidation  to 
C02  in  a  process  involving  the  transfer  of  two  protons  and  two 
electrons. 

The  system  that  fully  describes  the  CL  behavior  is  the  following: 
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1  Note  that  when  no  liquid  flux  to  the  CL  occurs,  this  intermittent  description 
excludes  the  effect,  but  not  the  presence,  of  liquid  accumulations  in  the  GDL 


ii  =  -7co)-exp(77/b1) 

Lref,l 


(12) 


C 

h  =  i*2-7^-7coexP  (V/b2)  (13) 

Eq.  (7)  describes  charge  conservation  in  the  CL,  while  Eqs.  (8) 
and  (9)  express  the  mass  balance  of  methanol  and  water,  respec¬ 
tively.  Due  to  the  strong  uncertainties  related  to  the  two-phase 
mass  transport  phenomena  in  the  CL  [43],  considering  for 
example  the  mixture  phase  transition  and  the  evolution  of  C02 
bubbles,  average  concentrations  and  diffusivities  are  considered  for 
simplicity.  In  particular,  the  average  diffusivities  of  methanol  and 
water  are  assumed  equal  to  each  other  and  proportional  to  the 
average  diffusivity  of  the  GDL: 


-'t.CHjOH 


=  D 


t,H20 


-r-^GDL.L  _  ^GDL.G  r  - 

AC  -sgdl-(-AC  -(1-Sgdl) 


■CF 


(14) 


CF  is  a  correction  factor,  that  takes  into  account  the  different 
porosity  and  void  fraction  of  the  CL. 

Eq.  (10)  is  the  Ohm’s  law  relating  the  proton  current  density 
to  the  gradient  of  overpotential;  it  is  necessary  to  take  into  ac¬ 
count  the  proton  transport  limitations  in  the  impedance  spec¬ 
trum.  Eq.  (11)  describes  the  variation  of  the  surface  coverage  by 
the  CO  adsorbed  intermediate:  the  first  reaction  step  produces 
CO  coverage,  while  the  second  reaction  step  consumes  CO 
coverage. 


3.3.  DMFC  anode  polarization  model 

During  anode  polarization  cathode  is  considered  a  reference 
electrode,  in  particular  a  dynamic  hydrogen  electrode  [44],  whose 
potential  is  nearly  constant  and  negligible  in  comparison  to  anode. 
Therefore  in  the  DMFC  anode  polarization  model  the  cathode 
overpotential  is  neglected  and  the  voltage  variations  are  mainly 
attributable  to  anode  and  membrane  overpotentials,  as  described 


M.  Zago,  A.  Casalegno  /  Journal  of  Power  Sources  248  (2014)  1181-1190 


1185 


in  Ref.  [45].  The  Eq.  (29)  in  Ref.  [34]  is  thus  replaced  by  the 
following: 

VCen  =  Va  +  Vohm  (15) 

Regarding  the  continuity  equations  at  cathode  side,  the 
following  modifications  have  been  introduced: 


J{S)  =  1.417  s- 2.120  s2  +  1.263  s3  (21) 

The  system  describing  CL  behavior  is  subject  to  the  following 
boundary  conditions: 

!Igdl  el  =  0  (22) 


•  air  is  replaced  by  hydrogen  and  therefore  Eqs.  (12)  and  (15)  in 
Ref.  [34]  are  substituted  with: 


!cell 


hc  a(vc-c£2)  _  i 

2  dx  2  F 


(16) 


-Dt 


a  Cr 


dx 


_  iwCH3OH 


(23) 


(24) 


•  CO2,  Eq.  (13)  in  Refs.  [34  ,  is  removed  because  the  methanol 
cross-over  flux  is  no  longer  oxidized; 

•  water  production  in  the  cathode  CL  is  not  considered  because  of 
the  hydrogen  feeding  and  Eq.  (14)  in  Ref.  [34]  becomes: 


hc  3(V-Ch20) 

2  ax 


JVH2° 

cross 


(17) 


Indeed  the  systems  characterizing  GDL,  Eqs.  (1)— (4),  and  CL 
behavior,  Eqs.  7—11,  have  been  introduced  in  the  above  described 
DMFC  anode  polarization  model.  The  GDL  equations  are  solved 
applying  the  appropriate  boundary  conditions  regarding  methanol 
concentration  and  saturation  at  the  channel-GDL  interface,  sch_GDL. 
The  complicated  two-phase  hydrodynamics  in  the  channel  makes 
the  liquid  saturation  at  this  interface  difficult  to  be  determined 
theoretically.  In  the  literature,  the  most  of  the  works  [46,47  ]  suppose 
its  value:  in  Ref.  [46]  it  is  equal  to  0.8,  while  in  Ref.  [47]  it  is  equal 
0.65,  despite  the  lower  hydrophobicity  of  the  considered  diffusion 
media.  Moreover,  the  above  cited  works  [46,47]  propose  a  ID  model 
along  the  thickness  of  the  DMFC  and  no  attempts  to  calculate  sch_GDL 
profile  along  channel  length  have  been  done.  In  Refs.  [48],  in  which  a 
complete  DMFC  model  is  presented,  the  saturation  at  channel-GDL 
interface  is  assumed  equal  to  the  channel  one:  however  this 
assumption  is  not  coherent  with  the  different  hydrophilic— hydro- 
phobic  properties  between  the  channel  and  the  diffusion  media.  In 
Refs.  [49],  instead,  an  empirical  approach  is  taken:  the  boundary 
condition  regarding  scii-gdl  profile  along  channel  length  is  expressed 
as  a  function  of  the  current  density.  According  to  Ref.  [49],  also  in  this 
work  an  empirical  approach  is  adopted.  However  the  saturation  is 
not  assumed  directly  proportional  to  the  current  density,  but  the 
following  dependence  on  saturation  in  the  channel  is  proposed: 

Sch-GDL  =  $1  'sch  +  S3  (18) 


where  sCh  is  the  saturation  in  the  channel  and  the  parameters  Si,  S2 
and  S3  are  calibrated  over  experimental  data.  Subsequently,  as 
widely  accepted  in  the  literature  [46,50],  the  liquid  saturation  at  the 
GDL-CL  interface  in  the  electrode  can  be  calculated  by  assuming 
that  the  capillary  pressure  remains  uniform  across  the  interface: 

PrlcDL-er  =  PclcDL-el+  (l9) 

According  to  Ref.  [35],  the  capillary  pressure  can  be  expressed 
as: 


Pc  =  <7H2o  ■  cosd c-(£/K)°  5  J(s)  (20) 

where  J(s)  is  the  Leverett  function  and  for  a  hydrophobic  medium  is 
given  by  the  following  relation: 


-Dt 


ac‘ 


0X 


=  N 


h2o 

GDL 


(25) 


-  ser^H30H|GDL_el  +  (1  -  sel)  ^-CH3OH | GDL_el 


(26) 


=  Sel'Qj 


2U  lGDL-el 


+  (1  ~~  sel)'CH20 


(27) 


This  is  a  two-point  boundary  value  problem  and  the  built  in  BVP 
function  of  MATLAB®  environment  has  been  used  for  the  numerical 
resolution  of  the  system  of  Eqs.  7—11. 

By  solving  the  entire  DMFC  anode  polarization  model  it  is 
therefore  possible  to  obtain  the  steady-state  profiles  along  channel 
length  and  both  GDL  and  CL  thickness  of  all  the  quantities  neces¬ 
sary  for  the  development  of  the  physically  based  DMFC  anode 
impedance  model,  reported  in  the  next  paragraph. 


3.4.  DMFC  anode  impedance  model 


After  obtaining  the  stationary  solution,  the  transient  perturba¬ 
tion  solution  of  the  DMFC  anode  can  be  obtained  in  the  following 
way.  Eqs.  (2)  and  (7)  are  linearized  about  a  steady-state  value  and 
subsequently  each  of  the  variables  z\  of  Eqs.  (1)— (4)  and  (7)— (11) 
are  perturbed  about  a  steady-state  value  with  a  sufficiently  low 
sinusoidal  disturbance: 


Zj(t)  =  z?  +  Re{Azrexp(/-w-t)}  (28) 

Neglecting  the  terms  with  products  of  the  disturbances  and 
subtracting  the  steady-state  equations,  it  is  possible  to  obtain  a 
system  of  linear  equations  for  the  complex  perturbation  amplitudes 
Azi  in  the  frequency  domain,  as  reported  in  Refs.  [20,21].  In 
particular  the  GDL  equations  (1)— (4)  become: 


9  AC? 


dx 


GDL,CH3OH  ’^GDL 


nG 

l,GDL,CH3OH 


■RT 


1-5 


.0  ' 

GDL  , 


rGDL,L,0 
LCH3OH 
‘  rGDL,L0 

LH20 


a  /-GDL.L 
nH2O,0  aLCH3OH 
iVni  ‘  rGDL,L,0 
LH20 


GDL 


-  Asgdl ■ 


L 

GDL,CH3OH 


nG 

uGDL,CH3OH 


ac, 


GDL,L  0 
CH3OH 


Kh.cHjOH-P-T/  dx 


(29) 
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9Asr 


Mh20'vh2o 


dX  tfHzO'COStfc-^GDL'^GDL)0'5  V  (SGDl) 


mH2O,0 

GDL  -Asr 


-  1  ■  ANH2° 

H%dl)  GDL 


SANC‘ 


+J-M'£  GDL 


_  j  ,  ,  c  f  Ac  /^GDL,L,0  ,  a /^GDL,L  c0  \ 

—  J-W-EGDL'  (^sGDL  '  cCH;iOH  +  a<-CH3OH'SGDLj 

•GDL,L,0 _ \  «"GDL,L  a^-GDL.L  o 

H  a'-CH3OH+aL'"»  ‘  “ 

%ch3oh-R-^ 


/Asgdl  •  GCH3qh  ^lch3oh+^lch3oh  SGDL 


aA/vH2C 

ULiiVGDL 

dx 

where: 


=  -J  ’ w  ’  fGDL  •  (  ^SGDL  •  CH  o 


■  f  AsGdl  ■  CuUnL’°  -  Asgdl  ■  C, 


.  ’sat  ' 
'GDL'lH20  , 


f{s gdl)  =  4DL- (T  -417  -4.24-sgdl  +  3.789-s£dl) 

While  the  CL  equations  7—11  take  the  following  form: 

9Ai  .  .  .  . 

—  =  Ai!  +Ai2+j-wCdl-/l?7 

0AN™3°h 
dx 


Ait  1  .  A  7=:t 

—  ^47p^J'w'Et'ACCH;iOH 


0AJVtH2°  M2  .  A  j 

- =  —  /  •  (0  ■  £t  •  ACo  n 

dx  2  F  J  1  H2° 


-D 


0AC 


ch3oh  =  anch3oh 


dx 


n  ^AC^j2q  _  ^fH20 

-Dt,H2°  -  A Nt 


dhri  A  i 
dx  <rt 


j-w-AyC0-,Tpt  = 


Aii  Ai2 


4-F  2F 
where  Aii  and  Ai2  are  equal  to: 

0 

'  ACch  oh 


a-  0>i  a  9,'i 

Aii  =  ^  -A V  +  ^r1 
a??l  0C 


a-  012 1 °  a  0i2 

Ai2  =  —  -An  h - ^ 

^1  0C 


CH3OH 

0 


9Tcol 


•^Tco 


h2o 


A7,t  di2  1°  a 

H2°  +  07^l  TC° 


(30) 


The  boundary  conditions  to  Eq.  29—32  and  Eq.  34- 
following: 

A^lcDL-el  =  ^EIS 


(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 

(41) 

(42) 
-40  are  as 

(43) 


Eq.  (43)  establishes  a  small  amplitude  perturbation  at  GDL-CL 
interface.  As  demonstrated  in  Refs.  [20],  the  system  of  Eq.  29—32 
and  Eq.  34—40  is  linear  and  homogenous  and  the  impedance 
does  not  depend  on  the  amplitude  of  this  boundary  condition. 
Eq.  (44)  means  null  proton  current  at  the  GDL-CL  interface.  Eq. 
(45)  expresses  zero  disturbance  of  methanol  and  water  concen¬ 
tration  at  channel-GDL  interface,  respectively.  This  implies  that 
the  concentration  in  the  channel  remains  unperturbed  by  the 
change  in  overpotential:  considering  the  high  liquid  accumula¬ 
tion  in  the  channel,  that  is  present  for  all  the  perturbation  fre¬ 
quencies,  this  is  a  reasonable  assumption.  Since  the  membrane  is 
permeable  to  both  methanol  and  water,  the  oscillation  of  the 
corresponding  flux  at  the  membrane  interface  cannot  be 
assumed  equal  to  zero.  As  reported  in  Refs.  [27],  water  cross-over 
flux  is  regulated  by  electro-osmotic  drag  and  liquid  diffusion  and 
convection,  while  the  methanol  cross-over  flux  is  only  governed 
by  electro-osmotic  drag  and  liquid  diffusion.  The  contributions  of 
liquid  diffusion  and  convection  are  influenced  by  the  corre¬ 
sponding  concentration  and  pressure  at  the  cathode  side,  that 
during  dynamic  operation  are  difficult  to  be  determined  with  a 
suitable  uncertainty.  Considering  also  the  considerable  accumu¬ 
lation  in  the  membrane,  Eq.  (46)  imposes  that  the  disturbance  of 
methanol  and  water  flux  is  equal  only  to  the  oscillation  of  drag 
component. 

Finally  the  impedance  of  the  anode  can  be  numerically  calcu¬ 
lated  as  the  ration  between  the  oscillating  voltage  and  current 
density  at  CL-membrane  interface: 


7  _  AtjI 

Ariel  m 


(47) 


The  system  of  Eq.  29—32  and  Eq.  34—40,  along  with  the 
boundary  conditions  Eq.  43—46,  is  solved  in  MATLAB®  environ¬ 
ment.  In  particular  the  FSOLVE  function  solves  the  CL  equations,  Eq. 
34—40,  by  iteratively  changing  the  value  of  the  oscillatory  fluxes  at 
electrode  inlet  till  the  boundary  conditions,  Eq.  (46),  are  satisfied. 
At  the  same  time  these  fluxes  are  imposed,  together  with  Eq.  (45), 
as  boundary  conditions  in  a  BVP  function,  that  solves  the  GDL 
equations,  Eq.  29—32.  In  this  way  all  the  boundary  conditions  re¬ 
ported  in  Eq.  43—46  are  simultaneously  satisfied  and  the  iterations 
are  performed  by  means  of  an  optimized  built  in  function, 
increasing  the  numerical  stability  of  the  problem. 

Furthermore,  considering  all  the  steady-state  profiles  obtained 
by  the  DMFC  anode  polarization  model,  it  is  possible  to  calculate 
the  local  value  of  DMFC  impedance  along  channel  length.  Thus  the 
total  impedance  is  equal  to  the  parallel  between  each  local 
impedance,  given  by: 


(48) 


where  N  is  the  number  of  integration  steps  in  which  the  channel 
length  has  been  discretized  during  the  numerical  resolution  of  the 
DMFC  anode  polarization  model. 


4.  Modeling  results  discussion 


A’lcDL-el  —  0 


ACCH30H|ch_GDL  =  0  ACH2o|ch  GDL  =  0 
ajvch3°h  I  =  AJVCH3°H  AJVH2°I 

lel-m  rn.drag  lel-m 


AJV 


H20 

m.drag 


^  4.1.  Steady-state  polarization 

The  increased  complexity  of  the  developed  model,  the  more 
(45)  detailed  description  of  the  DMFC  operation  and  the  necessity  to 
reproduce  both  steady-state  polarization  and  impedance  data 
make  the  calibration  procedure  more  difficult  than  in  the  previous 
models  [27,34],  However  the  model  reproduces  anode  polarization 
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Fig.  4.  Comparison  of  the  simulated  and  measured  anode  polarization. 

Table  1 

Assumed  and  calibrated  kinetic  and  mass  transport  parameters. 


DLGDL,CH30H 

10-5.4163-999.778/T  x  1Q4 

cm2  s-1 

[51] 

nG 

uGDL,CH3OH 

1.52  X  10~2 

cm2  s-1 

Calibrated 

Cref.l 

1  X  10~3 

mol  cm-3 

- 

Cref.2 

5.5  x  10~2 

mol  cm-3 

- 

ai 

0.74 

- 

Calibrated 

*2 

0.38 

- 

Calibrated 

h 

6.18  x  10-5  x  exp 
(6671 2/R  x  (1/353  -  1/T)) 

A  cm-3 

Calibrated 

1*2 

2.39  x  10~2  x  exp 
(6671 2/R  x  (1/353  -  1/T)) 

A  cm-3 

Calibrated 

PH20 

(-0.0028  x  (T  -  273)2 

gem-3 

[511 

-  0.1757  x  (T  -  273) 
+1003.8)  x  10~3 


MH20 

3.56  x  10-4 

Pa  s 

[511 

°H20 

6.25  x  10-2 

N  ITT1 

[51] 

Vt 

o 

X 

Q_1  cm-1 

Assumed 

^H,CH30H 

2.2  x  exp  (5200 
x  (1/T  -  1/298)) 

mol  J-1 

[51] 

Rgdl 

1  X  10-14 

m2 

Assumed 

6c 

130 

° 

Assumed 

F 

96485 

C  mol-1 

- 

R 

8314 

J  mor1  K1 

- 

/gdl 

3  x  10-2 

cm 

- 

ft 

o 

X 

cm 

- 

^cell 

22.1 

cm2 

- 

CF 

0.1 

- 

Calibrated 

S, 

0.75 

- 

Calibrated 

S2 

2 

- 

Calibrated 

$3,0.075  A/cm2 

0.085 

- 

Calibrated 

$3,0.15  A/cm2 

0.13 

- 

Calibrated 

Fig.  5.  Profile  of  sch  and  sch_CDL  along  channel  length  at  0.075  A  cm  2. 


Electrode  thickness  [%] 

Fig.  6.  Profile  of  current  density  along  channel  length  and  electrode  thickness  at 
0.075  A  cm-2. 


Table  2 

Assumed  parameters  for  the  simulation  of  anode  E1S. 


Cdi 

30 

F  cm-3 

A>t 

1  X  10-3 

mol  cm-3 

£t 

0.12 

- 

£gdl 

0.6 

- 

with  good  accuracy  (Fig.  4)  and  the  assumed  and  calibrated  pa¬ 
rameters  have  reasonable  values  (Table  1 ). 

The  high  number  of  calibration  parameters  reported  in  Table  1  is 
fundamental  to  simulate  also  the  anode  impedance  behavior:  in 


Fig.  7.  Comparison  of  the  simulated  and  measured  impedance  at  0.075  A  cm  2:  a) 
Nyquist  plot;  b)  Bode  plot. 


1188 


M.  Zago,  A.  Casalegno  /  Journal  of  Power  Sources  248  (2014)  1181-1190 


0.04 

0.03 


E 

N  0.01 


o  - 


-0.01 


Model  w/o  GDL 
'  Model  w  liquid  path 
Experimental 


0.02 


0.04 

Zreal  [ii] 


0.06 


0.08 


Fig.  8.  Simulated  anode  Nyquist  plot  at  0.075  A  cm  2  without  the  assumption  of 
intermittent  permeation  and  without  the  effect  of  GDL. 


fact  the  calibration  is  an  iterative  procedure  that  consists  in  mini¬ 
mization  of  the  residuals  between  model  estimation  and  experi¬ 
mental  measurements  of  both  anode  polarization  and  EIS.  A  limited 
number  of  calibration  parameters  would  be  necessary  to  reproduce 
only  the  polarization  curve. 

The  resulting  value  of  CF  is  lower  than  1  and  therefore  the  CL 
average  diffusivity  is  lower  than  the  diffusion  layer  one,  coherently 
with  reduced  values  of  porosity  and  void  fraction  in  the  CL. 
Regarding  the  calibration  parameters  of  saturation  profile  at 
channel-GDL  interface,  Eq.  (18),  the  values  of  Si  and  S2  are  quite 
high:  this  implies  a  very  strong  dependence  of  sCh-GDL  on  the 
saturation  in  the  channel.  Fig.  5  reports  a  comparison  between  the 
simulated  profile  of  sch  and  sCh-GDL  along  channel  length:  the 
saturation  at  the  channel-GDL  interface  is  lower  than  the  channel 
one,  coherently  with  the  higher  hydrophobicity  of  the  diffusion 
media. 

The  resolution  of  the  DMFC  anode  polarization  model  permits  to 
obtain  the  steady-state  profiles  of  all  the  physical  quantities 
necessary  for  simulation  of  impedance  behavior.  By  way  of  example 
Fig.  6  illustrates  current  density  profile  along  channel  length  and 
electrode  thickness:  its  value  is  zero  at  GDL  interface  and  reaches  its 
maximum  at  membrane  interface. 

4.2.  Electrochemical  impedance  spectroscopy 

4.2. ].  Impedance  simulations  at  0.075  A  cm~2 

All  the  kinetic  and  transport  parameters  necessary  for  the 
simulation  of  DMFC  anode  spectrum  have  been  previously  reported 


in  Table  1 ;  while  the  value  of  the  double  layer  capacitance,  active 
site  density  and  porosities  have  been  assumed,  Table  2. 

Fig.  7  a  reports  the  simulated  Nyquist  plot  at  0.075  A  cm-2.  The 
model  reproduces  experimental  observations  with  sufficient  ac¬ 
curacy  and  it  is  possible  to  notice  the  presence  of  the  linear  branch 
due  to  the  proton  transport  limitations  in  the  CL.  Moreover  the 
detailed  mass  transport  description  of  both  methanol  and  water 
fluxes  through  the  anode  GDL  provides  an  accurate  estimation  of 
the  mass  transport  impedance  and  the  total  resistance  is  close  to 
the  experimental  value.  Fig.  7b  illustrates  the  Bode  plot:  there  is  full 
agreement  between  model  and  experiment:  therefore  the  main 
involved  processes  are  correctly  described. 

It  is  very  important  to  figure  out  that  the  developed  model  re¬ 
produces  DMFC  operating  current  with  high  accuracy:  in  fact  the 
simulated  value  of  current  density,  equal  to  0.074  A  citT2,  is  almost 
identical  to  the  experimental  one.  This  is  a  very  important  feature 
because  it  would  be  possible  to  obtain  more  accurate  impedance 
simulations  by  slightly  changing  few  physical  parameters,  but  in 
this  case  the  DMFC  operating  current  would  be  rather  different 
from  the  experimental  one. 

Fig.  8  reports  impedance  simulations  performed  with  the 
assumption  that  GDL  is  always  flooded  with  fully  liquid  pathways: 
the  GDL  liquid  accumulation  terms  dampen  the  oscillations  of 
concentrations  for  all  the  perturbation  frequencies.  Model  pre¬ 
dictions  are  inconsistent  with  the  experimental  data  and  this  is  a 
further  confirmation  that  liquid  convection  through  the  GDL  is  an 
intermittent  phenomenon.  Furthermore  Fig.  8  illustrates  imped¬ 
ance  simulations  without  the  effect  of  the  GDL:  the  concentration 
disturbances  at  the  GDL-CL  interface  are  null  and  the  fluxes  are  not 
affected  by  the  GDL  accumulation  terms.  Even  at  low  current 
density  GDL  has  a  relevant  influence  on  impedance  behavior. 

Fig.  9  reports  a  comparison  between  model  results  with  and 
without  the  implementation  of  GDL  effects  along  channel  length.  It 
is  evident  that  towards  the  end  of  channel  the  mass  transport  losses 
due  to  the  presence  of  GDL  are  more  pronounced,  as  expected. 
Moreover,  also  the  magnitude  of  inductive  loop  increases  towards 
the  end  of  the  channel. 

Hereinafter  a  specific  analysis  is  performed  in  order  to  evaluate 
the  effects  of  the  boundary  conditions  regarding  the  oscillating 
fluxes  at  the  electrode-membrane  interface.  From  Fig.  10  it  is 
possible  to  figure  out  that  these  boundary  conditions  affect 
impedance  features  in  the  low-medium  frequencies  region.  In 
particular  ANCH3°H|  ,_m  =  0  implies  an  increase  of  the  total 
resistance  and  a  reduction  of  inductive  behavior,  while 
ANH2°|el  m  =  0  has  the  opposite  effect.  However  all  the  physical 
phenomena  are  interconnected:  switching  off  certain  physical 


Fig.  9.  Simulated  anode  Nyquist  plot  along  channel  length  at  0.075  A  cm  2  with  and 
without  GDL  effects. 


Fig.  10.  Simulated  anode  Nyquist  plot  at  0.075  A  cm  2  withANCH3°H|e|  m  =  0  and 
ANH2°|e,_m  =  0. 
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phenomena  may  affect  the  impedance  behavior  of  other  physical 
processes,  therefore  it  is  not  easy  to  find  out  the  origin  of  these 
different  impedance  features.  However  these  effects  seem  to  be 
related  to  a  variation  of  both  water  and  methanol  concentrations  in 
the  CL.  This  is  a  further  confirmation  of  the  complexity  of  DMFC 
anode  impedance  modeling  and  therefore  the  previous  model 
validation  on  three  different  typologies  of  measure  at  the  same 
time,  reported  in  Refs.  27],  was  fundamental  to  simulate  imped¬ 
ance  data  with  good  accuracy. 

4.2.2.  Impedance  simulations  at  0.15  A  cm~2 

Doubling  the  current  density  the  model  reproduces  experi¬ 
mental  observations  with  sufficient  accuracy,  Fig.  11:  the  inductive 
behavior  is  no  longer  present  and  mass  transport  limitations  seem 
to  appear  in  the  low  frequency  region.  However  the  model  slightly 
overestimates  the  total  resistance. 

In  Fig.  11,  the  simulations  without  the  effect  of  GDL  are 
considerable  different  from  the  experimental  observations.  As  ex¬ 
pected,  at  high  current  density  mass  transport  losses  are  more 
relevant  than  that  at  low  current  density,  Fig.  8. 

Once  again  the  model  predictions  with  the  effects  of  fully  liquid 
pathways  in  the  GDL  are  inconsistent  with  experimental  data.  In 
this  case  it  is  possible  to  clearly  distinguish  two  arches:  the  first  one 
characterizes  kinetic  losses,  while  the  second  one  is  peculiar  of 
mass  transport  limitations.  The  intermittent  liquid  convection  does 
not  dampen  the  oscillations  of  concentrations  at  GDL-CL  interface 
for  frequencies  lower  than  75  Hz,  leading  to  a  superimposition  of 
the  two  arches.  These  considerations  underline  the  difficulties  in 
the  interpretation  of  DMFC  anode  spectra:  in  fact  mass  transport 
limitations  do  not  manifest  themselves  as  a  second  arch,  even  if 
their  contribution  is  relevant. 

5.  Conclusions 

A  physically  based  model  of  DMFC  anode  impedance  has  been 
developed  and  validated  considering  two  different  operating  cur¬ 
rent  values.  The  model  includes  a  detailed  description  of  water  and 
methanol  transport  through  diffusion  layer:  the  former  is  regulated 
by  liquid  convection,  that  is  assumed  to  be  an  intermittent  phe¬ 
nomenon,  the  latter  is  governed  by  both  liquid  convection  and  gas 
and  liquid  diffusion.  Moreover  the  CL  has  a  finite  thickness  and 
methanol  oxidation  reaction  involves  CO  adsorbed  intermediate. 

The  developed  model  permits  to  quantify  DMFC  internal  losses 
and  investigate  the  origin  of  different  impedance  features.  The 
main  conclusions  regarding  DMFC  anode  impedance  behavior  are 
the  following: 


•  proton  transport  limitations  through  the  electrode  manifest 
themselves  as  a  45°  linear  branch  at  high  frequencies; 

•  the  intermittent  description  of  liquid  convection  through  GDL 
entails  that  liquid  accumulation  terms  dampen  the  oscillations 
of  concentrations  only  for  perturbation  frequencies  higher  than 
75  Hz.  This  assumption  implies  a  superimposition  between  the 
arches  peculiar  of  kinetic  and  mass  transport  losses  and  it  is  of 
fundamental  importance  to  correctly  simulate  anode  imped¬ 
ance  features; 

•  the  contribution  of  GDL  to  the  impedance  is  relevant  even  at  a 
rather  low  current  density,  such  as  0.075  A  cnT2,  and  increases 
along  channel  length,  as  expected.  Doubling  the  operating  cur¬ 
rent  density  the  effect  of  GDL  is  predominant  and  in  fact  without 
the  assumption  of  intermittent  permeation  it  is  possible  to 
clearly  distinguish  a  second  arch,  peculiar  of  mass  transport 
limitations: 

•  the  oscillations  of  methanol  and  water  fluxes  through  the 
membrane  cannot  be  neglected,  in  fact  they  affect  the  shape  of 
the  spectrum  in  the  low-medium  frequency  region:  this  in¬ 
creases  the  difficulties  in  the  validation  of  the  model  and  un¬ 
derlines  the  complexities  in  the  interpretation  of  DMFC  anode 
impedance.  A  specific  experimental  and  modeling  analysis  of 
the  effects  of  cross-over  fluxes  on  anode  impedance  behavior 
will  be  carried  out  in  future  work. 
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